In [1]:
import astropy.io.fits as fits
import numpy as np
from glob import glob
from matplotlib import pyplot as plt
import matplotlib.colors as colors
from astropy.visualization import simple_norm
import ccdproc as ccd
from astropy import visualization as viz
from astropy.visualization import (imshow_norm, MinMaxInterval, SqrtStretch, ZScaleInterval)
from astropy.coordinates import SkyCoord
from photutils.aperture import ApertureStats
from astropy.stats import SigmaClip
from photutils.aperture import aperture_photometry, CircularAperture, CircularAnnulus
from photutils.detection import DAOStarFinder
from astropy.stats import mad_std
from astropy.wcs import WCS
from astropy.table import Table
from astropy.io import ascii
from astropy import units as u
from datetime import datetime
from astropy.time import Time
In [2]:
# Sunset is 7:09 PM 
# By at observatory at 6:39 PM to take sky flats
# Ideally be done by 10 PM
In [3]:
data = ascii.read('landoltstars.txt')
In [4]:
star_pos = SkyCoord('16h41m41.24s', '36:27:35.5', unit=(u.hourangle, u.deg))
star_pos
Out[4]:
<SkyCoord (ICRS): (ra, dec) in deg
    (250.42183333, 36.45986111)>
In [5]:
data[0]
Out[5]:
Row index=0
StarDesignationRA(2000)Dec(2000)VB-VU-BV-RR-IV-I
str11str8str9float64float64float64float64float64float64
TPHE A00:30:09-46:31:2214.6510.7930.380.4350.4050.841
In [6]:
len(data)
Out[6]:
526
In [7]:
st = SkyCoord(data[0]['RA(2000)'], data[0]['Dec(2000)'], unit=(u.hourangle, u.deg))
st
Out[7]:
<SkyCoord (ICRS): (ra, dec) in deg
    (7.5375, -46.52277778)>
In [8]:
st.separation(star_pos).hour
Out[8]:
8.874273021712979
In [9]:
lands = []
for star in data:
    ra = star['RA(2000)']
    dec = star['Dec(2000)']
    coord = SkyCoord(ra, dec, unit=(u.hourangle, u.deg))
    #dist.append(coord.separation(star_pos))
    if star['V'] < 9.5 and coord.separation(star_pos).degree < 60:
        print(star)
        print('Separation:',coord.separation(star_pos).degree)
        lands.append(star)

StarDesignation RA(2000) Dec(2000)   V    B-V   U-B   V-R   R-I   V-I 
--------------- -------- --------- ----- ----- ----- ----- ----- -----
        109 231 17:45:20 -00:26:13 9.332 1.462 1.593 0.785 0.704 1.492
Separation: 39.74439851862938
StarDesignation RA(2000) Dec(2000)   V    B-V   U-B   V-R   R-I   V-I 
--------------- -------- --------- ----- ----- ----- ----- ----- -----
        111 773 19:37:17 +00:10:38 8.963 0.206 -0.21 0.119 0.144 0.262
Separation: 54.45338332240191
In [10]:
lands
Out[10]:
[<Row index=388>
 StarDesignation RA(2000) Dec(2000)    V      B-V     U-B     V-R     R-I     V-I  
      str11        str8      str9   float64 float64 float64 float64 float64 float64
 --------------- -------- --------- ------- ------- ------- ------- ------- -------
         109 231 17:45:20 -00:26:13   9.332   1.462   1.593   0.785   0.704   1.492,
 <Row index=431>
 StarDesignation RA(2000) Dec(2000)    V      B-V     U-B     V-R     R-I     V-I  
      str11        str8      str9   float64 float64 float64 float64 float64 float64
 --------------- -------- --------- ------- ------- ------- ------- ------- -------
         111 773 19:37:17 +00:10:38   8.963   0.206   -0.21   0.119   0.144   0.262]

Landolt Star Photometry¶

$V_{Landolt} = V_{Observed} + a (X-1) + b (X-1) (V_{Landolt}-I_{Landolt}) + ZP$

In [11]:
# For standard correction
# air mass
# airmass and color
# vega?
In [ ]:
In [ ]:
In [ ]:

Data Reduction¶

In [12]:
directory = '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/*'
glob(directory)
Out[12]:
['/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/bias',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/dark',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/flat_I',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/flat_V',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/Observation Plan (10_6).docx',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/ObservationPlan.pdf',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/vega']
In [13]:
files_bias = glob('/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/bias/*')
files_dark = glob('/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/dark/*')
files_flatI =  glob('/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/flat_I/*')
files_flatV =  glob('/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/flat_V/*')
#files_m2 =  glob('/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/*')
files_m2_I = glob('/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/*')
files_m2_V = glob('/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Vband/*')
files_vega_I =  glob('/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/vega/*_I_*')
files_vega_V =  glob('/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/vega/*_V_*')
In [14]:
def get_exptime(filename):
    return fits.getheader(filename)['EXPTIME']
def get_filter(filename):
    return fits.getheader(filename)['FILTER']
def get_obstime(filename):
    return Time(fits.getheader(filename)['DATE-OBS'])
def get_localtime(filename):
    return fits.getheader(filename)['NOTE'][11:]
In [15]:
names = ['Target', 'Exposure Time (s)', 'Filters', 'Number of Exposures', 'ObsTime (UTC)', 'Right Ascension', 'Declination']
filetype = ['Bias', 'Dark', 'I-Band Flat', 'V-Band Flat', 'Messier 2 I-Band', 'Messier 2 V-Band', 'Vega I-Band', 'Vega V-Band']
filelist = [files_bias, files_dark, files_flatI, files_flatV, files_m2_I, files_m2_V, files_vega_I, files_vega_V]
exptime = []
filters = []
numfiles = []
obstime = []
localtime = []
coord_m2 = SkyCoord('21h33m27s', '-00d49m23.7s', frame='fk5')
coord_vega = SkyCoord('18h36m56s', '+38d47m1s', frame='fk5')

coords_ra  = ['-', '-', '-', '-', '21h33m27s', '21h33m27s', '18h36m56s', '18h36m56s']
coords_dec = ['-', '-', '-', '-', '-00$^{\circ}$49\'23.7\"', '-00$^{\circ}$49\'23.7\"', '+38$^{\circ}$47\'1\"', '+38$^{\circ}$47\'1\"']

for lis in filelist:
    exptime.append(get_exptime(lis[0]))
    filters.append(get_filter(lis[0]))
    numfiles.append(len(lis))
    obstime.append(get_obstime(lis[0]))
    localtime.append(get_localtime(lis[0]))
In [16]:
Time(fits.getheader(files_m2_V[0])['DATE-OBS'])
Out[16]:
<Time object: scale='utc' format='isot' value=2022-10-07T01:51:05.000>
In [ ]:
In [17]:
fits.getheader(files_bias[0])['NOTE'][11:]
Out[17]:
'10/6/2022 at 18:57:46.000'
In [18]:
tblr = Table([filetype, exptime, filters, numfiles, obstime, coords_ra, coords_dec], names=names)
tblr.write('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive/obstable.dat', format='latex', overwrite=True)
In [19]:
pwd
Out[19]:
'/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive'
In [20]:
fits.getheader(files_bias[0])
Out[20]:
SIMPLE  =                    T /                                                
BITPIX  =                   16 /                                                
NAXIS   =                    2 /                                                
NAXIS1  =                  765 /                                                
NAXIS2  =                  510 /                                                
OBJECT  = '        '           /                                                
TELESCOP= 'Celestron CPC 12'''                                                  
INSTRUME= 'SBIG ST-402'                                                         
OBSERVER= 'Daniel_Emma_Alex'                                                    
DATE-OBS= '2022-10-06T22:57:46.000' / GMT START OF EXPOSURE [WIN]               
BZERO   = +3.276800000000E+004 /                                                
BSCALE  = +1.000000000000E+000 /                                                
EXPTIME = +4.000000000000E-002 / EXPOSURE IN SECONDS                            
CCD-TEMP= -5.232156845990E+000 / CCD TEMP IN DEGREES C                          
XPIXSZ  = +9.000000000000E+000 / PIXEL WIDTH IN MICRONS                         
YPIXSZ  = +9.000000000000E+000 / PIXEL HEIGHT IN MICRONS                        
XBINNING=                    1 / HORIZONTAL BINNING FACTOR                      
YBINNING=                    1 / VERTICAL BINNING FACTOR                        
XORGSUBF=                    0 / SUB_FRAME ORIGIN X POS                         
YORGSUBF=                    0 / SUB-FRAME ORIGIN Y POS                         
EGAIN   = +1.550000000000E+000 / ELECTRONS PER ADU                              
FOCALLEN= +3.048000000000E+003 / FOCAL LENGTH IN MM                             
APTDIA  = +3.048000000000E+002 / APERTURE DIAMETER MM                           
APTAREA = +6.840554060800E+004 / APERTURE AREA IN SQ-MM                         
CBLACK  =                 1044 / BLACK ADU FOR DISPLAY                          
CWHITE  =                 1173 / WHITE ADU FOR DISPLAY                          
PEDESTAL=                 -100 / ADD TO ADU FOR 0-BASE                          
DATAMAX =                65535 / SATURATION LEVEL                               
SBSTDVER= 'SBFITSEXT Version 1.0' / SBIG FITS EXTENSIONS VER                    
SWACQUIR= 'WinOPS Ver 5.61 Build 0-NT' / DATA ACQ SOFTWARE                      
SWCREATE= 'SBIG Win CCDOPS Version 5.61 Build 0-NT'                             
FILTER  = 'Photometric I'      / OPTICAL FILTER NAME                            
SNAPSHOT=                    1 / NUMBER IMAGES COADDED                          
DATE    = '2022-10-06'         / GMT DATE WHEN THIS FILE CREATED                
RESMODE =                    0 / RESOLUTION MODE                                
EXPSTATE= '125     '           / EXPOSURE STATE (HEX)                           
RESPONSE= +3.000000000000E+003 / CCD RESPONSE FACTOR                            
NOTE    = 'Local time:10/6/2022 at 18:57:46.000'
In [21]:
ccd_shape = fits.getdata(files_bias[0]).shape

Bias¶

In [22]:
bias = np.zeros((len(files_bias), ccd_shape[0], ccd_shape[1]))
for bb in range(len(files_bias)):
    #print(files_bias[bb])
    bias[bb] = fits.getdata(files_bias[bb])
In [23]:
master_bias = np.median(bias, axis=0)
plt.imshow(master_bias)
Out[23]:
<matplotlib.image.AxesImage at 0x7f2e06cee490>
In [24]:
master_bias.mean()
Out[24]:
1088.0907497116493
In [ ]:

Darks¶

In [25]:
exptime_darkI = fits.getheader(files_m2_V[0])['EXPTIME']
exptime_darkI
Out[25]:
20.0
In [26]:
exptime_darkV = fits.getheader(files_dark[0])['EXPTIME']
exptime_darkV
Out[26]:
20.0
In [27]:
dark = np.zeros((len(files_dark), ccd_shape[0], ccd_shape[1]))
for bb in range(len(files_dark)):
    #print(files_bias[bb])
    dark[bb] = fits.getdata(files_dark[bb])
In [28]:
master_dark = np.median(dark-master_bias, axis=0)
plt.imshow(master_dark)
np.mean(master_dark)
Out[28]:
9.130502370882994

Flats¶

In [29]:
exptime_flatI = fits.getheader(files_flatI[0])['EXPTIME']
In [30]:
flatI = np.zeros((len(files_flatI), ccd_shape[0], ccd_shape[1]))
for bb in range(len(files_flatI)):
    #print(files_bias[bb])
    flatI[bb] = fits.getdata(files_flatI[bb]) - master_bias - exptime_flatI/exptime_darkI * master_dark
In [31]:
master_flatI = np.median(flatI, axis=0)
plt.imshow(master_flatI)
np.mean(master_flatI)
Out[31]:
30056.210125195445
In [32]:
master_flatI_norm = master_flatI/np.mean(master_flatI)
plt.imshow(master_flatI_norm)
np.mean(master_flatI_norm)
Out[32]:
0.9999999999999997
In [33]:
exptime_flatV = fits.getheader(files_flatV[0])['EXPTIME']
In [34]:
flatV = np.zeros((len(files_flatV), ccd_shape[0], ccd_shape[1]))
for bb in range(len(files_flatV)):
    #print(files_bias[bb])
    flatV[bb] = fits.getdata(files_flatV[bb]) - master_bias - exptime_flatV/exptime_darkV * master_dark
In [35]:
master_flatV = np.median(flatV, axis=0)-master_bias
plt.imshow(master_flatV)
np.mean(master_flatV)
Out[35]:
25671.547239965388
In [36]:
master_flatV_norm = master_flatV/np.mean(master_flatV)
plt.imshow(master_flatV_norm)
np.mean(master_flatV_norm)
Out[36]:
1.0000000000000002

A Note about the flats:

These flats were taken as sky flats before astronomical twilight. As part of the process of aiming the telescope at the North star to center it, the CCD was removed from the back of the telescope. Afterwards, it was placed back on, but imperfectly so that the dust on the optics was not in the same place. In the calibrated science images to follow, the "donuts" from dust particles on the optics are visible because of this.

Calibration¶

$[(flat-bias)-\frac{t_f}{t_d}(dark-bias)]_{ave}[\frac{(object-bias)-\frac{t_o}{t_d}(dark-bias)}{(flat-bias)-\frac{t_f}{t_d}(dark-bias)}]_{pix}$

$[master\_flat-\frac{exptime\_flat}{exptime\_dark}(master\_dark)]_{ave}[\frac{(object-master\_bias) -\frac{exptime\_obj}{exptime\_dark}(master\_dark)}{(master\_flat)-\frac{exptime\_flat}{exptime\_dark}(master\_dark)}]_{pix}$

M2¶

In [37]:
files_m2_I
Out[37]:
['/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_001.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_002.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_003.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_004.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_001.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_002.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_003.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_004.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_005.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_006.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_007.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_008.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_009.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_010.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_011.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_012.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_013.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_014.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_015.FIT',
 '/mnt/c/Users/panda/Documents/ObsAstronomy/ObsData/LabFiveData/m2/Iband/m2_I20s_r2_016.FIT']
In [38]:
data_m2_I = np.zeros((len(files_m2_I), ccd_shape[0], ccd_shape[1]))

for ff in range(len(files_m2_I)):
    data_m2_I[ff] = (fits.getdata(files_m2_I[ff])-master_dark) / master_flatI_norm
In [39]:
pwd
Out[39]:
'/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive'
In [40]:
norm_ex = np.median(data_m2_I, axis=0)
imshow_norm(norm_ex, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
print('Very blurry because the telescope does not track well. Need to have something that accounts for the movement of the stars in every science frame.')
plt.tight_layout()
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive/plots/blurry.pdf')

Very blurry because the telescope does not track well. Need to have something that accounts for the movement of the stars in every science frame.
In [41]:
im1 = data_m2_I[0]
imshow_norm(im1, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
Out[41]:
(<matplotlib.image.AxesImage at 0x7f2df564f370>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2df56e06d0>)
In [42]:
bkg_sigma = mad_std(im1)  
daofind = DAOStarFinder(fwhm=4., threshold=3. * bkg_sigma)  
sources = daofind(im1)
for col in sources.colnames:  
    sources[col].info.format = '%.8g'  # for consistent table output
#print(sources)
In [43]:
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
apertures = CircularAperture(positions, r=5.)
In [44]:
fig = plt.figure(figsize=(8,9))
ax = imshow_norm(im1, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
apertures.plot(color='blue', lw=1.5, alpha=0.5)
#annulus_aperture.plot(color='blue', lw=1.5, alpha=0.1)
print('')










    
    































In [45]:



     


star = (677, 243)#(597, 34)



     



















In [46]:



     


#dist = np.sqrt((positions[:,0]-star[0])**2 + (positions[:,1]-star[1])**2)
#index = np.where(np.min(dist) == dist)
#pos = positions[index][0]
#star_pos.append(pos)
#ids.append(index)



     



















In [47]:



     


#pos



     



















In [48]:



     


from astropy.visualization import ImageNormalize



     























Image Combination¶


















In [49]:



     


#imshow_norm(im1, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
#shifted = []
data_m2_I = np.zeros((len(files_m2_I), ccd_shape[0], ccd_shape[1]))

for ff in range(len(files_m2_I)):
    im = (fits.getdata(files_m2_I[ff])-master_dark) / master_flatI_norm
    bkg_sigma = mad_std(im)  
    daofind = DAOStarFinder(fwhm=4., threshold=3. * bkg_sigma)  
    sources = daofind(im)
    positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
    apertures = CircularAperture(positions, r=5.) 
    dist = np.sqrt((positions[:,0]-star[0])**2 + (positions[:,1]-star[1])**2)
    index = np.where(np.min(dist) == dist)
    pos = positions[index][0]
    aper = CircularAperture(positions, r=5.) 
    #aper[index].plot(color='blue', lw=1.5, alpha=0.5)
    
    shi = np.roll(im, (int(star[1]-pos[1]), int(star[0]-pos[0])), axis=(0,1))
    
    #fig, axs = plt.subplots(1, 2, figsize=(9,9))
    #axs[0].imshow(im, norm=ImageNormalize(data=im, interval=ZScaleInterval(), stretch=viz.LinearStretch()), origin='lower', cmap='gray')
    #axs[0].scatter(pos[0], pos[1], marker='.', color='blue', alpha=0.5)
    #axs[0].scatter(star[0], star[1], marker='.', color='red', alpha=0.5)
    #axs[1].imshow(shi, norm=ImageNormalize(data=shi, interval=ZScaleInterval(), stretch=viz.LinearStretch()), origin='lower', cmap='gray')
    #axs[1].scatter(pos[0], pos[1], marker='.', color='blue', alpha=0.5)
    #axs[1].scatter(star[0], star[1], marker='.', color='red', alpha=0.5)
    
    data_m2_I[ff] = shi
    #print(im)
    
    #np.roll(im, pos-star)



     



















In [50]:



     


#star_V = (594, 290)
star_V = (555, 70)
#star_V = (355, 177)
#star_V = (600, 21)
#star_V = (599, 479)
#star_V = (682, 230)



     



















In [51]:



     


#imshow_norm(fits.getdata(files_m2_V[0]), origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
#shifted = []
data_m2_V = np.zeros((len(files_m2_V), ccd_shape[0], ccd_shape[1]))

for ff in range(len(files_m2_V)):
    im = (fits.getdata(files_m2_V[ff])-master_dark) / master_flatV_norm
    bkg_sigma = mad_std(im)  
    daofind = DAOStarFinder(fwhm=4., threshold=2.5 * bkg_sigma)  
    sources = daofind(im)
    positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
    apertures = CircularAperture(positions, r=5.) 
    dist = np.sqrt((positions[:,0]-star_V[0])**2 + (positions[:,1]-star_V[1])**2)
    index = np.where(np.min(dist) == dist)
    pos = positions[index][0]
    aper = CircularAperture(positions, r=5.) 
    #aper[index].plot(color='blue', lw=1.5, alpha=0.5)
    
    print(np.sqrt((pos[0]-star_V[0])**2 + (pos[1]-star_V[1])**2))
    shi = np.roll(im, (int(star_V[1]-pos[1]), int(star_V[0]-pos[0])), axis=(0,1))
    
    #fig, axs = plt.subplots(1, 2, figsize=(9,9))
    #axs[0].imshow(im, norm=ImageNormalize(data=im, interval=ZScaleInterval(), stretch=viz.LinearStretch()), origin='lower', cmap='gray')
    #axs[0].scatter(pos[0], pos[1], marker='.', color='blue', alpha=0.5)
    #axs[0].scatter(star_V[0], star_V[1], marker='.', color='red', alpha=0.5)
    #axs[1].imshow(shi, norm=ImageNormalize(data=shi, interval=ZScaleInterval(), stretch=viz.LinearStretch()), origin='lower', cmap='gray')
    #axs[1].scatter(pos[0], pos[1], marker='.', color='blue', alpha=0.5)
    #axs[1].scatter(star_V[0], star_V[1], marker='.', color='red', alpha=0.5)
    
    if np.sqrt((pos[0]-star_V[0])**2 + (pos[1]-star_V[1])**2) > 30.:
        print('Selected star too far away')
    else: 
        data_m2_V[ff] = shi
    #print(im)
    
    #np.roll(im, pos-star)



     






















    
    






26.53738398585733
21.83304490504877
21.953036079431303
46.23193614817681
Selected star too far away
12.23206935111877
10.70735201945291
9.572564077172094
9.091323261693075
8.902041990812997
43.06073330763313
Selected star too far away
1.3256166289841604
5.831531623950198
11.455285389002357
5.18858963701437
5.318678335272187
5.538186868950662
14.841000883137074
14.090627098947493
11.151455903061647
9.690312857984582
























In [52]:



     


image_combined_I = np.median(data_m2_I, axis=0)
fig = plt.figure(figsize=(8,9))
ax = imshow_norm(image_combined_I, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')



     






















    
    































In [53]:



     


fig = plt.figure(figsize=(8,9))
ax = imshow_norm(data_m2_I[0], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')



     






















    
    































In [54]:



     


image_combined_V = np.median(data_m2_V, axis=0)
fig = plt.figure(figsize=(8,9))
ax = imshow_norm(image_combined_V, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')



     






















    
    































In [55]:



     


fig = plt.figure(figsize=(8,9))
ax = imshow_norm(data_m2_V[0], origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')



     






















    
    































In [56]:



     


bkg_sigma = mad_std(image_combined_I)  
daofind = DAOStarFinder(fwhm=4., threshold=3. * bkg_sigma)  
sources = daofind(image_combined_I)
for col in sources.colnames:  
    sources[col].info.format = '%.8g'  # for consistent table output
    
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
apertures = CircularAperture(positions, r=5.)  
annuli = CircularAnnulus(positions, 5., 8.)

fig = plt.figure(figsize=(8,9))
ax = imshow_norm(image_combined_I, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
apertures.plot(color='blue', lw=1.5, alpha=0.2)
annuli.plot(color='blue', lw=1.5, alpha=0.2)

print('')



     






















    
    















    
    































In [57]:



     


aperstats = ApertureStats(image_combined_I, annuli)
bkg_mean = aperstats.mean
phot_table_I = aperture_photometry(image_combined_I, apertures)
aperture_area = apertures.area_overlap(image_combined_I)
total_bkg = bkg_mean * aperture_area
phot_bkgsub = phot_table_I['aperture_sum'] - total_bkg
phot_table_I['total_bkg'] = total_bkg
phot_table_I['aperture_sum_bkgsub'] = phot_bkgsub



     



















In [58]:



     


bkg_sigma = mad_std(image_combined_V)  
daofind = DAOStarFinder(fwhm=4., threshold=3. * bkg_sigma)  
sources = daofind(image_combined_V)
for col in sources.colnames:  
    sources[col].info.format = '%.8g'  # for consistent table output
    
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
apertures = CircularAperture(positions, r=5.)  
annuli = CircularAnnulus(positions, 5., 8.)

fig = plt.figure(figsize=(8,9))
ax = imshow_norm(image_combined_V, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
apertures.plot(color='blue', lw=1.5, alpha=0.2)
annuli.plot(color='blue', lw=1.5, alpha=0.2)

print('')



     






















    
    















    
    































In [59]:



     


aperstats = ApertureStats(image_combined_V, annuli)
bkg_mean = aperstats.mean
phot_table_V = aperture_photometry(image_combined_V, apertures)
aperture_area = apertures.area_overlap(image_combined_V)
total_bkg = bkg_mean * aperture_area
phot_bkgsub = phot_table_V['aperture_sum'] - total_bkg
phot_table_V['total_bkg'] = total_bkg
phot_table_V['aperture_sum_bkgsub'] = phot_bkgsub



     



















In [60]:



     


#color = phot_table_V['aperture_sum_bkgsub'] - phot_table_I['aperture_sum_bkgsub']
#magnitude = phot_table_V['aperture_sum_bkgsub']



     



















In [61]:



     


x_VIshift = phot_table_V[4]['xcenter'] - phot_table_I[4]['xcenter']
y_VIshift = phot_table_V[4]['ycenter'] - phot_table_I[4]['ycenter']



     



















In [62]:



     


x_VIshift



     






















    
    
Out[62]:







$5.1888917 \; \mathrm{pix}$























In [63]:



     


y_VIshift



     






















    
    
Out[63]:







$-12.35472 \; \mathrm{pix}$























In [64]:



     


int(star_V[1]-pos[1])



     






















    
    
Out[64]:








-8
























In [65]:



     


r = np.roll(image_combined_V, (-int(x_VIshift.value), -int(y_VIshift.value)), axis=(1,0))
plt.imshow(r)



     






















    
    
Out[65]:








<matplotlib.image.AxesImage at 0x7f2df42fa7f0>










    
    



































Photometry¶


















In [66]:



     


bkg_sigma = mad_std(r)  
daofind = DAOStarFinder(fwhm=4., threshold=3. * bkg_sigma)  
sources = daofind(r)
for col in sources.colnames:  
    sources[col].info.format = '%.8g'  # for consistent table output
    
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
apertures = CircularAperture(positions, r=5.)  
annuli = CircularAnnulus(positions, 5., 8.)

fig = plt.figure(figsize=(8,9))
ax = imshow_norm(r, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
apertures.plot(color='blue', lw=1.5, alpha=0.2)
annuli.plot(color='blue', lw=1.5, alpha=0.2)

print('')



     






















    
    















    
    































In [67]:



     


aperstats = ApertureStats(image_combined_V, annuli)
bkg_mean = aperstats.mean
phot_table_V = aperture_photometry(image_combined_V, apertures)
aperture_area = apertures.area_overlap(image_combined_V)
total_bkg = bkg_mean * aperture_area
phot_bkgsub = phot_table_V['aperture_sum'] - total_bkg
phot_table_V['total_bkg'] = total_bkg
phot_table_V['aperture_sum_bkgsub'] = phot_bkgsub



     



















In [68]:



     


nearest_star = np.zeros(len(phot_table_V))
for st in range(len(phot_table_V)):
    xV = phot_table_V[st]['xcenter']
    yV = phot_table_V[st]['ycenter']
    nearest = 999
    
    for ar in range(len(phot_table_I)):
        xI = phot_table_I[ar]['xcenter']
        yI = phot_table_I[ar]['ycenter']
        dist = np.sqrt((xV-xI)**2+(yV-yI)**2).value
        
        if dist < nearest:
            nearest_star[st] = int(ar)
            nearest = dist



     



















In [69]:



     


nearest_star



     






















    
    
Out[69]:








array([  0.,   1.,   2.,   3.,   4.,   5.,   6.,   7.,   8.,  15.,   9.,
        15.,  10.,  11.,  12.,  13.,  15.,  16.,  17.,  20.,  21.,  22.,
        25.,  24.,  25.,  26.,  27.,  28.,  30.,  29.,  31.,  32.,  30.,
        31.,  33.,  34.,  35.,  36.,  37.,  38.,  39.,  42.,  41.,  44.,
        43.,  45.,  46.,  49.,  48.,  50.,  51.,  52.,  54.,  55.,  59.,
        57.,  60.,  58.,  61.,  62.,  63.,  64.,  65.,  66.,  74.,  69.,
        67.,  68.,  71.,  70.,  72.,  74.,  73.,  75.,  76.,  82.,  78.,
        77.,  81.,  80.,  83.,  84.,  85.,  86.,  87., 105.,  89.,  90.,
        91.,  92.,  94.,  93.,  95.,  96.,  97.,  96.,  99.,  98., 103.,
       102., 104., 105., 109., 106., 111., 107., 108., 110., 112., 114.,
       113., 116., 115., 118., 120., 117., 121., 119., 122., 123., 124.,
       125., 115., 126., 126., 127., 127., 128., 129., 130., 131., 133.,
       132., 136., 134., 135., 138., 137., 139., 141., 140., 142., 151.,
       143., 144., 146., 147., 148., 144., 149., 150., 151., 152., 153.,
       154., 155., 157., 156., 158., 159., 160., 162., 161., 163., 164.,
       165., 167., 166., 168., 169., 170., 171., 172., 173., 175., 174.,
       177., 179., 180., 168., 182., 184., 183., 185., 186., 189., 188.,
       190., 191., 192., 194., 195., 196., 197., 198., 200., 201., 202.,
       203., 204., 205., 206., 207., 208.])
























In [70]:



     


#color = phot_table_V['aperture_sum_bkgsub'] - phot_table_I['aperture_sum_bkgsub']
#magnitude = phot_table_V['aperture_sum_bkgsub']



     



















In [71]:



     


nearest_star.size



     






















    
    
Out[71]:








204
























In [72]:



     


mag = []
color = []
for ss in range(nearest_star.size):
    xs = phot_table_V[ss]['xcenter'] - phot_table_I[int(nearest_star[ss])]['xcenter']
    ys = phot_table_V[ss]['ycenter'] - phot_table_I[int(nearest_star[ss])]['ycenter']
    dist = np.sqrt(xs**2 + ys**2).value
    #print(dist)
    if dist < 1.0:
        mag.append(phot_table_V[ss]['aperture_sum_bkgsub'])
        color.append(phot_table_V[ss]['aperture_sum_bkgsub'] - phot_table_I[int(nearest_star[ss])]['aperture_sum_bkgsub'])



     



















In [73]:



     


plt.scatter(color, mag)



     






















    
    
Out[73]:








<matplotlib.collections.PathCollection at 0x7f2df4123e50>










    
    































In [74]:



     


image_combined_I.shape



     






















    
    
Out[74]:








(510, 765)
























In [75]:



     


fig = plt.figure(figsize=(8,9))
ax = imshow_norm(image_combined_I, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
apertures.plot()
CircularAperture((295, 310), r=125.).plot(color='blue', alpha=0.5)



     






















    
    
Out[75]:








(<matplotlib.patches.Circle at 0x7f2df406aa30>,)










    
    































In [76]:



     


mask = np.ones(image_combined_I.shape)
mask_bool = np.ones(image_combined_I.shape)
center = (295, 310)
radius = 125.
for x in range(mask.shape[1]):
    for y in range(mask.shape[0]):
        dist = np.sqrt((center[0] - x)**2+(center[1] - y)**2)
        if dist <= radius:
            mask[y,x] = 0#np.nan
            mask_bool[y,x] = np.nan
        dist2 = np.sqrt((mask.shape[1]/2 - x)**2+(mask.shape[0]/2 - y)**2)
        if dist2 >= 450:
            mask[y,x] = 0#np.nan
            mask_bool[y,x] = np.nan
        #print(dist)

plt.imshow(mask)



     






















    
    
Out[76]:








<matplotlib.image.AxesImage at 0x7f2dde8fc100>










    
    































In [77]:



     


fig = plt.figure(figsize=(8,9))
ax = imshow_norm(r, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
CircularAperture((295, 310), r=125.).plot(color='blue', alpha=0.5)
CircularAperture((mask.shape[1]/2, mask.shape[0]/2), r=400.).plot(color='blue', alpha=0.5)



     






















    
    
Out[77]:








(<matplotlib.patches.Circle at 0x7f2dde8e3310>,)










    
    































In [78]:



     


im_final_V = r*mask
im_final_I = image_combined_I*mask
im_final_I



     






















    
    
Out[78]:








array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       ...,
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.]])
























In [79]:



     


bkg_sigma = mad_std(im_final_V)  
daofind = DAOStarFinder(fwhm=4., threshold=3.0 * bkg_sigma)  
sources = daofind(im_final_V)
for col in sources.colnames:  
    sources[col].info.format = '%.8g'  # for consistent table output
    
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
apertures = CircularAperture(positions, r=5.)  
annuli = CircularAnnulus(positions, 5., 8.)

fig = plt.figure(figsize=(8,9))
ax = imshow_norm(im_final_V*mask_bool, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
apertures.plot(color='blue', lw=1.5, alpha=0.2)
annuli.plot(color='blue', lw=1.5, alpha=0.2)

aperstats = ApertureStats(im_final_V, annuli)
bkg_mean = aperstats.mean
phot_table_V = aperture_photometry(im_final_V, apertures)
aperture_area = apertures.area_overlap(im_final_V)
total_bkg = bkg_mean * aperture_area
phot_bkgsub = phot_table_V['aperture_sum'] - total_bkg
phot_table_V['total_bkg'] = total_bkg
phot_table_V['aperture_sum_bkgsub'] = phot_bkgsub

print('')
plt.tight_layout()
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive/plots/photometry.pdf', bbox_inches='tight')



     






















    
    















    
    































In [80]:



     


bkg_sigma = mad_std(im_final_I)  
daofind = DAOStarFinder(fwhm=4., threshold=3.0 * bkg_sigma)  
sources = daofind(im_final_I)
for col in sources.colnames:  
    sources[col].info.format = '%.8g'  # for consistent table output
    
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
apertures = CircularAperture(positions, r=5.)  
annuli = CircularAnnulus(positions, 5., 8.)

fig = plt.figure(figsize=(8,9))
ax = imshow_norm(im_final_I*mask_bool, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
apertures.plot(color='blue', lw=1.5, alpha=0.2)
annuli.plot(color='blue', lw=1.5, alpha=0.2)

aperstats = ApertureStats(im_final_I, annuli)
bkg_mean = aperstats.mean
phot_table_I = aperture_photometry(im_final_I, apertures)
aperture_area = apertures.area_overlap(im_final_I)
total_bkg = bkg_mean * aperture_area
phot_bkgsub = phot_table_I['aperture_sum'] - total_bkg
phot_table_I['total_bkg'] = total_bkg
phot_table_I['aperture_sum_bkgsub'] = phot_bkgsub

print('')



     






















    
    















    
    































In [81]:



     


nearest_star = np.zeros(len(phot_table_V))
for st in range(len(phot_table_V)):
    xV = phot_table_V[st]['xcenter']
    yV = phot_table_V[st]['ycenter']
    nearest = 999
    
    for ar in range(len(phot_table_I)):
        xI = phot_table_I[ar]['xcenter']
        yI = phot_table_I[ar]['ycenter']
        dist = np.sqrt((xV-xI)**2+(yV-yI)**2).value
        
        if dist < nearest:
            nearest_star[st] = int(ar)
            nearest = dist



     



















In [82]:



     


mag = []
color = []
for ss in range(nearest_star.size):
    xs = phot_table_V[ss]['xcenter'] - phot_table_I[int(nearest_star[ss])]['xcenter']
    ys = phot_table_V[ss]['ycenter'] - phot_table_I[int(nearest_star[ss])]['ycenter']
    dist = np.sqrt(xs**2 + ys**2).value
    #print(dist)
    if dist < 1.0:
        mag.append(phot_table_V[ss]['aperture_sum_bkgsub'])
        color.append(phot_table_V[ss]['aperture_sum_bkgsub'] - phot_table_I[int(nearest_star[ss])]['aperture_sum_bkgsub'])



     



















In [83]:



     


plt.scatter(color, mag)
plt.gca().invert_yaxis()



     






















    
    































In [84]:



     


fig = plt.figure(figsize=(8,9))
ax = imshow_norm(r, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
# Final V-band
plt.tight_layout()
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive/plots/messier2V.pdf', bbox_inches='tight')



     






















    
    































In [85]:



     


fig = plt.figure(figsize=(8,9))
ax = imshow_norm(image_combined_I, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')
# Fianl I-Band
plt.tight_layout()
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive/plots/messier2I.pdf', bbox_inches='tight')



     






















    
    



































Vega¶






















$[master\_flat-\frac{exptime\_flat}{exptime\_dark}(master\_dark)]_{ave}[\frac{(object-master\_bias)
-\frac{exptime\_obj}{exptime\_dark}(master\_dark)}{(master\_flat)-\frac{exptime\_flat}{exptime\_dark}(master\_dark)}]_{pix}$























Image Combination¶


















In [86]:



     


exptime_vegaI = fits.getheader(files_vega_I[0])['EXPTIME']



     



















In [87]:



     


exptime_darkI



     






















    
    
Out[87]:








20.0
























In [88]:



     


data_vega_I = np.zeros((len(files_vega_I), ccd_shape[0], ccd_shape[1]))

for ff in range(len(files_vega_I)):
    data_vega_I[ff] = (fits.getdata(files_vega_I[ff])-(exptime_vegaI/exptime_darkI)*master_dark) / master_flatI_norm



     



















In [89]:



     


vega_final_I = np.median(data_vega_I, axis=0)
imshow_norm(vega_final_I, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')



     






















    
    
Out[89]:








(<matplotlib.image.AxesImage at 0x7f2df43e4940>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2ddeac5eb0>)










    
    































In [90]:



     


exptime_vegaV = fits.getheader(files_vega_V[0])['EXPTIME']



     



















In [91]:



     


data_vega_V = np.zeros((len(files_vega_V), ccd_shape[0], ccd_shape[1]))

for ff in range(len(files_vega_V)):
    data_vega_V[ff] = (fits.getdata(files_vega_V[ff])-(exptime_vegaV/exptime_darkV)*master_dark) / master_flatV_norm



     



















In [92]:



     


vega_final_V = np.median(data_vega_V, axis=0)
imshow_norm(vega_final_V, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')



     






















    
    
Out[92]:








(<matplotlib.image.AxesImage at 0x7f2df42a4d00>,
 <astropy.visualization.mpl_normalize.ImageNormalize at 0x7f2df4396fd0>)










    
    



































Note: Vega is definitely over exposed























Photometry¶


















In [93]:



     


bkg_sigma = mad_std(vega_final_I)  
daofind = DAOStarFinder(fwhm=15., threshold=100.0 * bkg_sigma)  
sources = daofind(vega_final_I)
for col in sources.colnames:  
    sources[col].info.format = '%.8g'  # for consistent table output
    
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
apertures = CircularAperture(positions, r=60.)  
annuli = CircularAnnulus(positions, 5., 8.)

ap_I = CircularAperture((498, 197), r=60.)
an_I = CircularAnnulus((498, 197), 70., 120.)

fig = plt.figure(figsize=(8,9))
ax = imshow_norm(vega_final_I, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')

ap_I.plot(color='red', lw=1.5, alpha=0.2)
an_I.plot(color='red', lw=1.5, alpha=0.2)

#apertures.plot(color='blue', lw=1.5, alpha=0.2)
#annuli.plot(color='blue', lw=1.5, alpha=0.2)

print('')
plt.tight_layout()
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive/plots/vegaI.pdf', bbox_inches='tight')



     






















    
    















    
    































In [94]:



     


aperstats = ApertureStats(vega_final_I, an_I)
bkg_mean = aperstats.mean
phot_table_vegaI = aperture_photometry(vega_final_I, ap_I)
aperture_area = ap_I.area_overlap(vega_final_I)
total_bkg = bkg_mean * aperture_area
phot_bkgsub = phot_table_vegaI['aperture_sum'] - total_bkg
phot_table_vegaI['total_bkg'] = total_bkg
phot_table_vegaI['aperture_sum_bkgsub'] = phot_bkgsub

phot_table_vegaI



     






















    
    
Out[94]:







QTable length=1



idxcenterycenteraperture_sumtotal_bkgaperture_sum_bkgsub


pixpix


int64float64float64float64float64float64


1498.0197.059047290.5342532412600461.47632765846446829.05792558


























In [ ]:



     


 



     



















In [95]:



     


bkg_sigma = mad_std(vega_final_V)  
daofind = DAOStarFinder(fwhm=16., threshold=100.0 * bkg_sigma)  
sources = daofind(vega_final_V)
for col in sources.colnames:  
    sources[col].info.format = '%.8g'  # for consistent table output
    
positions = np.transpose((sources['xcentroid'], sources['ycentroid']))  
apertures = CircularAperture(positions, r=80.)  
annuli = CircularAnnulus(positions, 5., 8.)

ap_V = CircularAperture((469, 185), r=75.)
an_V = CircularAnnulus((469, 185), 100., 150.)

fig = plt.figure(figsize=(8,9))
ax = imshow_norm(vega_final_V, origin='lower', interval=ZScaleInterval(), stretch=viz.LinearStretch(), cmap='gray')

ap_V.plot(color='red', lw=1.5, alpha=0.2)
an_V.plot(color='red', lw=1.5, alpha=0.2)

#apertures.plot(color='blue', lw=1.5, alpha=0.2)
#annuli.plot(color='blue', lw=1.5, alpha=0.2)

print('')
plt.tight_layout()
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive/plots/vegaV.pdf', bbox_inches='tight')



     






















    
    















    
    































In [96]:



     


aperstats = ApertureStats(vega_final_V, an_V)
bkg_mean = aperstats.mean
phot_table_vegaV = aperture_photometry(vega_final_V, ap_V)
aperture_area = ap_V.area_overlap(vega_final_V)
total_bkg = bkg_mean * aperture_area
phot_bkgsub = phot_table_vegaV['aperture_sum'] - total_bkg
phot_table_vegaV['total_bkg'] = total_bkg
phot_table_vegaV['aperture_sum_bkgsub'] = phot_bkgsub

phot_table_vegaV



     






















    
    
Out[96]:







QTable length=1



idxcenterycenteraperture_sumtotal_bkgaperture_sum_bkgsub


pixpix


int64float64float64float64float64float64


1469.0185.0115364750.2482222219942595.75966265895422154.48855956


























In [97]:



     


phot_table_I



     






















    
    
Out[97]:







QTable length=82



idxcenterycenteraperture_sumtotal_bkgaperture_sum_bkgsub


pixpix


int64float64float64float64float64float64


113.5540624491241811.84660995481611561986.2596530730653826.4551919042358159.804461168824


2751.42504616221451.851170779937263262041.9759183886353734.784572164738307.191346223895


3595.105232628340433.10370087945881100503.6093102833192804.31963971457699.289670568818


420.4199239654261158.7927925184922394673.8273802052692551.692449671782122.1349305334734


5256.862278139844365.62829102346362100167.175921722693532.44994334186634.725978380811


6337.987900399974477.3539535715546898702.9891333534393535.271153988435167.717979364999


7549.81648725534482.2000357425378298768.8839769422892853.51134682925915.372630113081


8230.9234083973787285.27939424051158101669.1238743027294218.681882451137450.441991851592


9333.474586226284890.52576834219148101319.2145367883393921.228386626727397.986150161611


10188.58014055544743121.733561987702496492.310133189193359.354259875563132.9558733135345


11512.7921102699793123.5917214558294797293.9326109540792840.212824662224453.719786291855


12349.57846761791666130.8166465939787113854.2375308289596632.5482751621217221.68925566683


13129.77302367955207131.243645256748199218.1019996338795290.84426289763927.2577367362683


14122.77009753307367135.189988216222497495.912116424895454.490094742322041.4220216824906


15278.1797807816903136.9606640195692496182.0940462297393668.099786375722513.9942598540074


..................


68204.03202978693398428.9174337003223699248.4169759978597790.328597025521458.0883789723302


69400.54870368739836428.943209378363798833.6354217654795432.283201607543401.3522201579326


70194.48150990383192433.9323301665131124338.46285111597100554.7607036125823783.702147503383


71195.09763830666935442.9130690818934799777.950489144100878.42210998737-1100.4716208433674


72137.49988682444481451.0578691384308101620.3573150909796701.389555810794918.9677592801745


73196.35965946954394457.9901744254795104313.7860582777696346.212431275857967.5736270019115


74280.9762331964428459.04722896464483125933.56636823025100638.4465978058925295.119770424368


75343.9128415509698460.0612440892582599254.243788084695425.64136571783828.602422366792


76133.54062267889867461.3818727820196102670.4946463463396426.590974357696243.903671988635


77622.8511800857611474.369425225922598102.7264624399294515.393909973143587.332552466789


78445.0018572895247483.227263620490599748.4186423931894133.628355144175614.7902872490085


79327.9762114437748488.519061703318102915.016255579595292.449203705667622.567051873833


80290.11404339062364489.8215848313834102748.7518605241895774.42100429216974.330856232074


81592.551732093303491.3400594929506107048.315400201795480.1053367809411568.21006342076


82186.57758077001452499.57231970824216112927.4800970787396145.6584477462716781.821649332458






























Magnitude Calculations¶
m$_1 - $ m$_{ref} = -2.5 log_{10}(\frac{I_1}{I_{ref}})$



















In [98]:



     


exptime_m2I = exptime_darkI
exptime_m2V = exptime_darkV
gain = fits.getheader(files_m2_I[0])['EGAIN'] # electrons per adu



     



















In [99]:



     


intensity_m2V = phot_table_V['aperture_sum_bkgsub']*gain/exptime_m2V
intensity_m2I = phot_table_I['aperture_sum_bkgsub']*gain/exptime_m2I

intensity_vegaV = phot_table_vegaV['aperture_sum_bkgsub']*gain/exptime_vegaV
intensity_vegaI = phot_table_vegaI['aperture_sum_bkgsub']*gain/exptime_vegaI



     



















In [100]:



     


mag_m2V = -2.5*np.log10(intensity_m2V/intensity_vegaV)
mag_m2I = -2.5*np.log10(intensity_m2I/intensity_vegaI)



     






















    
    






/tmp/ipykernel_5060/3524271707.py:2: RuntimeWarning: invalid value encountered in log10
  mag_m2I = -2.5*np.log10(intensity_m2I/intensity_vegaI)
























In [101]:



     


phot_table_V['magnitude'] = mag_m2V
phot_table_I['magnitude'] = mag_m2I



     



















In [102]:



     


nearest_star = np.zeros(len(phot_table_V))
for st in range(len(phot_table_V)):
    xV = phot_table_V[st]['xcenter']
    yV = phot_table_V[st]['ycenter']
    nearest = 999
    
    for ar in range(len(phot_table_I)):
        xI = phot_table_I[ar]['xcenter']
        yI = phot_table_I[ar]['ycenter']
        dist = np.sqrt((xV-xI)**2+(yV-yI)**2).value
        
        if dist < nearest:
            nearest_star[st] = int(ar)
            nearest = dist



     



















In [103]:



     


mag_V = []
mag_I = []
color = []
for ss in range(nearest_star.size):
    xs = phot_table_V[ss]['xcenter'] - phot_table_I[int(nearest_star[ss])]['xcenter']
    ys = phot_table_V[ss]['ycenter'] - phot_table_I[int(nearest_star[ss])]['ycenter']
    dist = np.sqrt(xs**2 + ys**2).value
    #print(dist)
    if dist < 1.0:
        mag_V.append(phot_table_V[ss]['magnitude'])
        mag_I.append(phot_table_I[int(nearest_star[ss])]['magnitude'])
        color.append(phot_table_V[ss]['magnitude'] - phot_table_I[int(nearest_star[ss])]['magnitude'])



     



















In [104]:



     


fig, axs = plt.subplots(1, 2, figsize=(9,5))

axs[0].scatter(color, mag_I)
axs[0].invert_yaxis()
axs[0].set_xlabel('Color (V-I)')
axs[0].set_ylabel('I-band Magnitude')

axs[1].scatter(color, mag_V)
axs[1].invert_yaxis()
axs[1].set_xlabel('Color (V-I)')
axs[1].set_ylabel('V-band Magnitude')

plt.tight_layout()
plt.savefig('/mnt/c/Users/panda/Documents/HomeworkBack/ObsTech/LabFive/plots/colormag.pdf')



     






















    
    































In [ ]:



     


 



     



















In [ ]: